A model for complex aftershock sequences 
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The decay rate of aftershocks is commonly very well described by the modified Omori law, 
n(t) ex t~ p , where n(t) is the number of aftershocks per unit time, t is the time after the main 
shock, and p is a constant in the range 0.9 < p < 1.5, and usually close to 1. But there are 
also more complex aftershock sequences for which the Omori law can be considered only as a first 
approximation. One of these complex aftershock sequences took place in the Eastern Pyrenees 
on February 18, 1996, and was described in detail by Correig et al. [1997]. In this paper, we 
propose a new model inspired by dynamic fiber-bundle models to interpret this type of complex 
aftershock sequences with sudden increases in the rate of aftershock production not directly related 
to the magnitude of the aftershocks (as in the epidemic-type aftershock sequences). The model is a 
simple, discrete, stochastic fracture model where the elements (asperities or barriers) break because 
of static fatigue, transfer stress according to a local load-sharing rule and then are regenerated. We 
find a very good agreement between the model and the Eastern Pyrenees aftershock sequence and 
we propose that the key mechanism for explaining aftershocks, apart from a time-dependent rock 
strength, is the presence of dynamic stress fluctuations which constantly reset the initial conditions 
for the next aftershock in the sequence. 
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I. INTRODUCTION 

Omori discovered scaling in earthquakes in the fre- 
quency distribution of aftershocks over one hundred years 
ago when he proposed a formula to represent the decay 
of aftershock activity with time [Omori, 1894]. Now, one 
hundred years later, it remains as one of the few well es- 
tablished empirical laws in Seismology. As noted by Utsu 
[1995], 'any theory for the origin of aftershocks must ex- 
plain this law, which is unique for its power law depen- 
dence on time'. The Omori law (as modified by Utsu, 
1961), 



n(t) = Kt~ 



(1) 



says that the number of aftershocks n(t), measured at 
time t after the time of the main shock, declines fol- 
lowing a power law with an exponent p around one 
(0.9 < p < 1.5, with a median of about 1.1, [Utsu, 1995]), 
and K being a proportionality constant. To avoid diver- 
gence at t = 0, the Omori law is usually written in the 
form 



n(t) = K(t + c)~ 



(2) 



where c is an additional 'small' positive constant with 
dimensions of time (between 0.01 and 1 days, with a me- 
dian of 0.3 days, [Utsu, 1995]) . The power law, scale-free 
behavior is maintained for t ^> c, with a transition to 
n(t) = const for t < c that accounts for incompleteness 
in the detection of low-magnitude aftershocks during the 



few hours after the mainshock (e.g. [Utsu, 1995, fig 2; 
Gross and Kisslinger 1994]. From Eq. (|^), the cumula- 
tive number of aftershocks, N(t), occurred until time t 
after the main shock, defined as J Q n(s)ds, is 



N(t) = 



_ ( K\n(t/c+l) 



if p = 1 



\K{< 



i-p _ 



(i + c )l-P}/(p-l) if p? 4 1 



(3) 



The Omori law has also been verified in laboratory- 
scale experiments in brittle rock deformation by measur- 
ing acoustic emission [Scholz, 1968a, b; Lockner and By- 
erlee, 1977; Hirata, 1987; Sammonds et al, 1992; Lock- 
ner, 1993] and in mine-induced seismicity [Talabi, 1997], 
which represents an intermediate scale between lab ex- 
periments and natural aftershocks. The fulfilment of the 
Omori law at the microscale (10 -3 — 10 _1 m), mesoscale 
(1 — 10 1 m) and macroscale (10 2 — 10 4 m), suggests that a 
universal process is behind the inelastic strain responsi- 
ble for acoustic emission in the laboratory [Hirata, 1987], 
induced microseismicity in mines [Gibowicz, 1997] and af- 
tershock sequences in active tectonic faults [Utsu, 1995; 
Gross and Kisslinger 1994]. But, what is this mecha- 
nism? 

Static fatigue, also known as stress-, creep-, or delayed 
fracture is the basic way of time-dependent failure under 
a constant load of a broad variety of materials, including 
textile fibers [Coleman, 1957], fiber composites [Phoenix, 
1977], wood [Garcimartin et al., 1997], microcrystals 
[Pauchard and Meunier, 1993], gels [Bonn et al., 1998], 
policrystalline ceramics [Jacobs and Chen, 1994], metals 
[Schleinkofer et al, 1996], silicate glasses [Charles, 1958], 
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minerals [Scholz, 1968a; Barnett and Kerrich, 1980], and 
rocks [Atkinson, 1984]. In all these cases, the signature 
of static fatigue is the observation of a failure strength 
that is a function of the load history of the material. For 
brittle materials, and from the point of view of fracture 
mechanics, time-dependent strength is commonly asso- 
ciated with kinetic fracture, i.e., with the propagation 
of cracks under a crack tip stress intensity factor below 
the modulus of cohesion of the material [Kostrov et at, 
1969]. This propagation is stable and quasi-static, and 
is referred to as subcritical crack growth, where 'quasi- 
static' means at velocities much less than the sonic ve- 
locity of the medium [Das and Scholz , 1991]. The pres- 
ence of a chemically active fluid environment saturating 
the pore and crack space enhances this subcritical crack 
growth, a mechanism known as stress corrosion [Charles, 
1958; Wiederhorn, 1967]. There is ample evidence that 
geological materials under brittle conditions owe their 
time-dependent strength to the mechanism of subcriti- 
cal crack growth assisted by stress corrosion [Atkinson, 
1984; Atkinson and Meredith, 1987]. 

Benioff [1951] presented the first detailed theory offer- 
ing an explanation of the causes and characteristics of 
aftershock sequences in terms of identifiable mechanical 
properties. According to his theory, aftershocks occur 
when there is a time- dependent recovery of stress follow- 
ing the main shock. The stress recovery was ascribed by 
Benioff to a static fatigue of the rocks in the immediate 
area of the fault. 

Since this seminal paper, many laboratory and numer- 
ical experiments have confirmed the hypothesis that af- 
tershocks are a process of relaxing stress concentrations 
produced by the dynamic rupture of the main shock, and 
that they are, therefore, an intrinsic time-dependent rhe- 
ological effect. In this context, Scholz [1968b] formulated 
the first time-dependent strength model of aftershocks. 
He suggested that a time-dependent strength of the rocks 
in the area of the main shock could be the cause of the 
aftershock sequences and invoked static fatigue due to lo- 
cal overloads to stresses much higher than their long-term 
strength as the main mechanism of aftershocks. Based 
on Scholz's [1968a] laboratory experiments on static fa- 
tigue of quartz, Das and Scholz [1981] formulated a gen- 
eral model of aftershocks using elastic fracture mechanics 
and the concept of subcritical crack growth. They showed 
that this model is consistent with the decay rate of after- 
shocks as expressed by the Omori law, and that it is able 
to reproduce many other characteristics of real aftershock 
sequences. More recent works and papers that stress the 
role of time-dependent strength in aftershock dynamics 
are: Yamashita and Knopoff [1987], who assume that 
stress corrosion is the physical mechanism for the delayed 
fracture in aftershocks; Marcellini [1995, 1997], advocat- 
ing static fatigue, together with stress inhomogeneities, 
as the cause of Omori-law aftershock sequences; and Lee 
[1999] and Lee and Sornette [1999], who constructed a 
fuse network model of aftershocks incorporating a time 
dependent strength compatible with the mechanism of 



subcritical crack growth. All these models of aftershocks 
also obey the Omori law. 

As mentioned above, the Omori-law decay rate of af- 
tershocks following a main shock is an almost universal 
characteristic of seismicity (as compared to the more ir- 
regular patterns of premonitory activity as foreshocks or 
quiescence). But despite this universality, many real af- 
tershock sequences display anomalies in the decay rate 
that depart from the simple Omori-law behavior. Among 
these anomalies we can cite [Utsu, 1995]: (i) cases in 
which seismic activity following the main shock cannot 
be represented by a simple power law due to the mixing 
of different series of activity [Gross and Kisslinger 1994]; 
and (ii) cases where aftershocks decay, as a whole, ac- 
cording to the Omori law, but depart temporarily from 
the formula due to abrupt changes in activity (accelera- 
tions and/or quiescence). 

In this paper we are interested in aftershock series that 
do not rigorously follow the Omori law and, in particu- 
lar, in this second type of anomalies where sudden ac- 
celerations in the rate of aftershock activity are not di- 
rectly linked to aftershocks of larger magnitude. This last 
case is the so-called epidemic-type aftershock sequence, 
ETAS, where each aftershock has its own sequence of 
aftershocks [Ogata, 1988], and can be thought of as a 
fractal version of the simple Omori relaxation formula. 
There are, however, some aftershock sequences where the 
changes in decay rate are independent of the magnitude 
of the aftershocks that provoke these changes in activity, 
and that cannot be ascribed to the ETAS model. One 
of these aftershock sequences took place in the eastern 
Pyrenees on February 18, 1996 [Correig et ai, 1997] and 
in this paper we propose a simple stochastic model a la 
dynamic fiber bundle model as a framework to interpret 
this class of aftershock sequences. 



II. THE DATA 

On February 18, 1996, a local magnitude Ml — 5.2 
earthquake occurred in the Eastern Pyrenees, with epi- 
central location N42deg47.71', E2deg32.30' and focal 
depth of 8 km [Rigo et ai, 1997; Pauchet et at, 1999]. 

The series of aftershocks that followed this event was 
recorded at the three-component continuous broadband 
seismic station at the Tunel del Cadi, located at about 80 
km SW of the epicentral area. Altogether, the series con- 
sists of 337 events (complete for a threshold magnitude 
of 1.9), spanning 1846 hours (77 days) from the time of 
the main shock, and with magnitudes ranging from 1.9 
to 3.8. Figure [l]a shows the cumulative series of after- 
shocks, along with the magnitude of the events. The 
sudden change in slope at about 300 hours is not due to 
incompleteness of the series, and from the point of view 
of the magnitude of the aftershocks, there is no specific 
characteristic, nor any relevant event, that justifies this 
sudden change in the event rate. Because of this diffcr- 
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ent behavior, we will restrict our attention to the series 
defined by the first 300 hours, with a total of 308 events, 
as displayed in Figure 0b. 

The most striking feature of this series is the change in 
concavity of the cumulative curve not correlated to any 
significant event (as would be the case from the point 
of view of an ETAS model), suggesting an increase in 
the rate of aftershocks production apparently not related 
to any relaxation process. If we try to fit by the maxi- 
mum likelihood method the aftershock data in Fig. |l|b to 
Eq. (3), we immediately appreciate that a unique fit to 
the whole range (0-300 hours) is graphically worse than 
two independent fits to the ranges 0-100 hours and 140- 
300 hours. Unfortunately, the parameter estimation for 
the 140-300 hours is not unique, in the sense that more 
than one set of parameter values return the same value of 
the likelihood function after maximization. After fixing 
the values of K and c (to the values K = 17 ± 2 and 
c = 0.2 ±0.1 obtained for the fit to the 0-100 hours inter- 
val), the 140-300 hours fit returns a unique value for the 
parameter p = 0.654 ± 0.005, which is statistically differ- 
ent to the value p = 0.75 ± 0.04 obtained for the 0-100 
hours interval. The error in both estimates is one stan- 
dard deviation assuming a normal distribution for the 
variance. These values of p are abnormally low (0.75 for 
the 0-100 hours interval, and 0.65 for the 140-300 hours 
interval). The fit and the value of p do not improve if the 
magnitude threshold is increased [Correig et al., 1997]. 

The interpretation of the Omori law as a relaxation 
process suggests a way of separating the aftershocks in 
the series into two classes: class A for the events that 
follow a relaxation law and class B for those events that 
do not. The criterion to assign the events to classes A 
or B is the following: if the interval of time Ati between 
events i and i — 1 is strictly larger than the interval of 
time Ati-i between events i — 1 and i — 2, the event i 
belongs to class A; otherwise it belongs to class B. Events 
belonging to class A are termed leading aftershocks, and 
those belonging to class B, cascades. Figure |l|c shows the 
aftershock sequence classified as leading events (solid cir- 
cles) and cascades (dots) . Note that a cascade is initiated 
by a leading aftershock and that this leading aftershock 
has no significative different magnitude. 

The fit of Eq. (3) to the series formed by the leading 
aftershocks is shown in Fig. ||a. The fit has improved 
significantly from the initial fit to the whole sequence, 
and the value obtained for the exponent is now p = 0.94, 
much more in agreement with the standard values for 
worldwide aftershock sequences. Figure ||b depicts the 
series of cascades, in which the first term of each cascade 
is a leading aftershock. Two important features are read- 
ily visible from the figure: (i) the cascades are in general 
well approximated by straight lines; and (ii) their corre- 
sponding slopes decrease with time. A plot of the slope 
s of the cascades against time t (Figure ||) shows the re- 
markable fact that there exists a power-law relationship 
of the form s oc t~ v between them, with v w 0.71. 

The properties summarized in Figs. [j] through || for the 



aftershock sequence of the Eastern Pyrenees can be de- 
scribed at first order with the modified Omori law, Eqs. 
(2) and (3). But at second order there are important 
non-random fluctuations about this law (represented by 
the cascades) that can not be fitted in detail with, nor 
accounted for, the Omori law and its relaxation origin. 
In the next Section we will construct a model to account 
for this second order deviations from the Omori law, and 
for the Omori law itself, of course. 

We want to stress here that the characteristics of the 
series of aftershocks from the February 18, 1996, Pyre- 
nees mainshock are by no means 'exceptional'. On the 
contrary, they seem to be a rather general feature of 
aftershock series. The authors are currently analyzing 
various aftershock sequences (Greece, Kobe, Landers, 
Northridge) and have found a behavior very similar to 
that of the Pyrenees aftershock sequence. Results will be 
reported elsewhere. 



III. DYNAMIC FIBER-BUNDLE MODELS 

Fiber-bundle models (FBMs) are simple discrete 
stochastic fracture models amenable to either close ana- 
lytical or fast numerical solution. These models arose in 
intimate connection with the strength of bundles of tex- 
tile fibers [Daniels, 1945; Coleman, 1957]. Since Daniels' 
and Coleman's seminal works there has been a long tradi- 
tion in the use of these simple models to analyze failure of 
heterogeneous materials [Vazquez- Prada et al., 1999 and 
references therein]. 

The dynamic version of the FBM simulates the failure 
of materials because of static fatigue or delayed rupture. 
In this version, one considers: (i) a discrete set of N 
elements located on the sites of a <i-dimensional lattice; 
(ii) a probability distribution for the nominal lifetimes of 
individual elements; and (Hi) a load-transfer rule which 
determines how the load carried by a failed element is to 
be distributed among the surviving elements in the set. 

As stated in (ii), the nominal lifetimes, tj, of the in- 
dividual elements supporting an initial stress Oi, equal 
for all j, are taken from a probability distribution of the 
type 

n 3 =l-e- k ^ j = 1,2,..., TV. (4) 

where nj are random numbers (0 < nj < 1) and k(o~i) 
is the so called hazard rate or breaking rule. The most 
accepted hazard rate is of the form 

k(a i ) = V o(^) f ', (5) 

Uq and <7o represent a hazard rate of reference and a 
stress of reference respectively, and p is an exponent in 
the range 2 < p < 50. This function has been used 
to fit experimental results of time-to-failure on various 
materials [Coleman, 1957; Phoenix, 1977 ]. Besides, 
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Phoenix and Tierney [1983] derived it from a kinetic the- 
ory of thermally activated atomic bond rupture [Zhurkov, 
1965], and showed that in many circumstances it is a bet- 
ter approximation than the exponential breaking rule, 
k(u) = aexp(/3cr), also used in modeling time dependent 
fracture [Coleman, 1957]. 

Equation (||) has the same form as Charles power- 
law to describe stress corrosion induced subcritical crack 
growth in geological materials [Atkinson, 1984]: 



v = v exp(-H/RT)Kf, 



(6) 



where v is the crack velocity, H the activation energy, 
R the gas constant, T the absolute temperature, Kj the 
stress intensity factor for mode I fracture, and vq and n 
are constants. Sometimes, n is known as the stress cor- 
rosion index. Nominal values at room temperature and 
in wet rock are [Atkinson and Meredith, 1987]: 15-40 for 
quartz and quartz rocks; 10-30 for calcite rocks; 30-70 
for granitic rocks; and 25-50 for gabbro and basalt. If we 
assume constant temperature, Eq. (^) can be simplified 
to 



v = AKf, 



(7) 



which is identical to Eq. (||) if we substitute a by K and 
identify the breaking rate expressed by Eq. (|5|) with the 
crack opening velocity expressed by Eq. (^) . 

In dynamic FBMs, once elements begin to fail because 
of fatigue and their stress is transferred according to the 
assumed rule, the stresses among the surviving elements 
are no longer equal and the stress history of an individual 
element becomes complicated by the successive step-like 
transfers coming from failing elements. The effect of the 
increase in stress for a particular element j is the reduc- 
tion of its lifetime from the initially assigned tj to a Tj 
defined by 



U 



dt. 



(8) 



Notice that in the case of independent elements (i.e., 
no stress transfer), (Tj(t) — Oi Vj and hence the bundle 
would break as a succession of individual failures at the 
times tj assigned at the beginning. When, on the con- 
trary, stress redistribution between elements is assumed, 
the temporal series of individual failures actually occurs 
at the times Tj dictated by Eq. (Sh. 

Note, that in these models, the total stress acting on 
the system is conserved until the failure of the last ele- 
ment; and for the reasons explained above, the last ele- 
ment, i.e., that with the longest Tj, does not coincide, in 
general, with that with the longest tj. A detailed expla- 
nation of how to perform a Monte Carlo simulation using 
this model can be found in [Newman et ai, 1995] where 
the reader is referred to for details. 

The dynamic FBM can also be calculated using a prob- 
abilistic approach. This approach was introduced by us 
in Gomez et al, [1998]. From this perspective, one starts 



with TV elements loaded with an initial common stress 
equal to o~i. The mean time interval, S, for one element 
to break by fatigue is 



1 



Eli 



(9) 



Supposing that v$ = 1 and o-$ = 1 in Eq. (Q), we have 

1 



5 = 



p ' 
'3=1 °i 



(10) 



In the first step, o~j — at Vj and hence S = j^pr . This 
will change with time because of stress transfers. At any 
instant of the process of breaking, 5 of Eq. (|l(]) represents 
the mean time for the next individual failure. The identi- 
fication of which element breaks after one S is calculated 
by deciding that the probability that precisely element k 
be the affected one is given by 



Pk 



(11) 



The reader will note that, from the probabilistic per- 
spective, one does not consider weak elements and strong 
elements: here all elements are equal but, in general, with 
a different o~j . The succession of individual breakings pro- 
ceeds by chance with the probabilities dictated by Eq. 
( p"l| ) until the total collapse of the system. In Gomez et 
ai, [1998] it was shown that the probabilistic approach 
represents a way of partially smoothing the fluctuations 
inherent to these stochastic models of fracture. 



IV. THE MODEL 

The model used in this paper to describe aftershock 
sequences is based on the same physical mechanism as 
the model recently introduced by Lee [1999] and Lee and 
Sornette [1999], although the ingredients of the cellular 
automaton, the way of running the model and the type 
of results are different. Our model is inspired by dy- 
namic FBMs, but with several substantive differences. 
As stated above, in dynamic FBMs N uniformly loaded 
elements break by fatigue one after one, in a total stress 
conserving process until the last individual failure. The 
sequence always occurs in a finite-time accelerated pro- 
cess. 

Trying to describe a sequence of aftershocks, which is 
a clearly decelerated process of relaxation, we will mod- 
ify three things: (i) stress is lost in the breaking process 
in two ways, one is by stress transfers out of the system 
through the borders, and the other way is by consider- 
ing a dissipative effect during each transfer. When an 
element fails bearing a stress a, the fraction (1 — 7r)er 
is removed from the system. Thus, the constant n will 
represent the degree of conservation. 

(ii) In FBMs, when an element fails, it transfers its 
load and then remains inactive. In this model, on the 
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contrary, when an element fails, it transfers its load (ex- 
cept the dissipated fraction) to its nearest neighbors but 
then it is automatically regenerated and able to receive 
stress again and actively participate in the time evolu- 
tion of the set. After its regeneration, an element has 
zero load before receiving any load. This assumption is 
justified because time intervals between individual fail- 
ures are much longer than the time of actual breaking 
of an asperity, and hence in an interevent interval a just 
failed asperity has enough time to reheal. 

The third difference (Hi) with FBMs is related to the 
initial stress distribution in the system. Whilst in FBMs 
it is usually assumed that the initial load per element 
is a constant o~i, here, trying to simulate the disordered 
state existing in the fault after a main shock, we will 
take the initial Oj from a uniform probability distribu- 
tion (0 < <Tj < 1; 1 < j < N). In the actual running of 
the model we will adhere to the probabilistic approach 
explained above for the FBM, and so we will deal with 
the <5s as defined in Eq. (nw and the probabilities of Eq. 
&■ 

Besides the three general modifications (z), (n), and 
(iii) introduced with respect to FBMs, in order to rein- 
force the appearance of sudden accelerations, which con- 
stitute the genuine phenomenology of the complex after- 
shocks sequences, we will add two more rules for running 
our model. 

When a breaking is going to occur in a context such 
that all the elements have a < 1, we say that this is a 
normal event, the 5 is calculated using Eq. (1C) and the 
broken element is pointed out by using Eq. (11). On the 
contrary if a breaking is going to occur with one or sev- 
eral elements with a > 1, we call it an avalanche step. 
For each avalanche step, a) the 5 is calculated using, as 
usual, Eq. ( |l0| ) but the element that breaks is that one 
whose <jj surpasses 1. If there are several elements with 
a > 1, the element that fails is that with the maximum 
value of o-j. b) The avalanche ends when all the Uj of the 
set become lower than 1. During an avalanche, which in 
general involves several <5s, all the elements that have sur- 
passed at any step of the avalanche the condition a > 1 , 
remain inactive with a = until the end of the avalanche. 

These two additional rules a) and b) are introduced 
in order to increase the local stress accumulations. The 
high stress concentrations occurring in avalanche events 
lead to very short <5s and in very short 5s it is reasonable 
to assume that there is not enough time for the healing 
process. 

As a resume of the last paragraphs, we will recall that 
in the process of evolution of the system there are normal 
events and avalanche events. The former (normal) refers 
to the failure of one element when no element in the sys- 
tem has a > 1. The latter (avalanche steps) corresponds 
to the failure of one element with a > 1. Eq. ([h]) is al- 
ways used for the calculation of the time intervals. The 
5s of the avalanche steps are much shorter because of the 
large stress concentrations induced by rules a) and b), 
and the magnitude of the exponent p. With these rules, 



it is obvious that the avalanches become extinct with 
time because as the total stress in the system declines, 
it is more difficult to locally accumulate load to surpass 
unity. The model of Lee [1999] and Lee and Sornette 
[1999] cannot describe complex aftershock sequences be- 
cause in their model the avalanches are instantaneous in 
time, so that the decay rate would have singularities at 
the time of an avalanche. On the contrary, our avalanches 
are not instantaneous since they are formed by several 
steps with their corresponding <5s, that is, they occur in 
a finite time interval. 

We have performed our simulations in a two dimen- 
sional square lattice with 50a;50 sites with p = 30 and 
7r = 0.7. Each site represents the emplacement of one 
element or asperity. The load-transfer rule assumed is 
a local load-sharing (LLS) rule: a failing element trans- 
fers its load to its nearest four neighbors located in the 
North, South, East and West. If the element is located 
at the borders of the square lattice this isotropic transfer 
provokes the corresponding stress leakage. 

In FBMs, in which broken elements remain inactive, an 
LLS rule can lead to an uncertain situation when a failing 
element lacks active neighbors able to receive the stress 
transfer. In this model, though, this does not occur be- 
cause any element, at any time, is active to receive stress. 
The only possible exceptional situation could come from 
the application of rule b) in avalanches. Then, the load 
of the failing element is transferred only to the existing 
active nearest neighbors. In this extremely rare case (we 
have not met such a situation in our simulations) of hav- 
ing the four nearest neighbors already broken, the load 
is removed from the lattice. This is just one possibility 
among various choices. Lomnitz-Adler et al. [1992] ex- 
plored in detail three of these possibilities for a cellular 
automaton with rules similar to those of the LLS static 
fiber-bundle model, and the reader is referred to this pa- 
per for details. In our model deciding among one of the 
three alternative scenarios is not important due to the 
extreme rarity of such events. 

The running of the program proceeds as follows. At 
t = all the elements on the lattice are initialized by 
loading them with a random initial stress o~j taken from 
a uniform probability distribution < <jj < 1. Then, we 
calculate the 5 corresponding to the first failure by using 
Eq. ( |To| ) . The choice of the element that actually breaks 
is done by using Eq. ([ll]) , and a random number between 
and 1 to materialize the choice. The chosen element 
fails and 7r times its stress is transferred to its nearest 
neighbors and the (1 — it) fraction disappears. Now, we 
analyze the distribution of stress in the board; if all the 
<7j are lower than one, the process of calculating the next 
failure is identical to the first. If, on the contrary, one 
(or several) <jj > 1 then we calculate the corresponding 
value of 5 from Eq. jl0| ) and the failure is assigned to the 
element with the biggest <jj > 1. During the period of 
an avalanche, rule b) is applied to favor the stress con- 
centration. Thus, the series of breakings and transfers, 
involving normal events or avalanches, proceeds until a 
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prescribed minimum value for the total stress on the sys- 
tem is reached, otherwise there will be an infinite number 
of aftershocks. 

In our model, due to the dissipation, the total stress 
in the system, S = ^2j&j, systematically decreases. If 
the value p — 1 were considered, then from Eq. (10), 



the successive 5s would necessarily be longer and longer. 
But p is bigger than one, and this is the reason why one 
can have a step down in S and find a shorter value of 
5. This is the key point to understand our model and 
other models based on subcritical crack growth. In the 
general trend of S reduction and hence temporal decel- 
eration, the stress transfers in the system provoke local 
inhomogeneities in a, and due to the high values of p, this 
leads to temporal accelerations. These accelerations are 
embedded in the general trend of dynamical relaxation. 



V. RESULTS AND CONCLUSIONS 

We have carried out numerical simulations which show 
the fulfillment of Omori's law and reproduce the features 
already commented on, that is, a cumulative plot with 
sudden variations in the number of events (accelerations). 
We show here the results for a two-dimensional system 
of 50 x 50 elements located on a square lattice, with p 
equal to 30 and a conservation level of tt — 0.7. Other 
simulations have also been performed varying the size of 
the system, the value of p and the conservation level, tt. 
We have also explored a variant of rule b) in which, in 
avalanches, all the elements with a > 1 break simulta- 
neously in the same 5. We have found that the results 
are indeed very close to those exposed here with equal 
qualitative behavior. Nevertheless, it should be noted 
that although the results are robust over a large range 
of parameters, different characteristics arise for extreme 
values of p and tt. 

Figure |] shows the rate of aftershocks dN/dt as a func- 
tion of time. Time is represented in dimensionless units 
and it is the sum of the successive 5s. The straight line 
has a slope of —1 for comparison. Thus, the 1/t decay 
is confirmed and is in full agreement with Omori's law 
for real aftershock sequences. The power law depicted 
is very robust over a wide range of the parameters that 
characterize the model. The most critical parameter is 
the conservation level since for values of tt close to unity, 
the system does not dissipate enough to avoid its com- 
plete failure. Besides, for large dissipation tt <C 1, the 
power law extends only to a few decades, and the num- 
ber of decades decreases as tt decreases. Anyhow, in all 
cases the exponent of the power law decay is very close 
to unity. The major vertical spikes of Fig. ^ correspond 
to avalanche-type events that disappear with time. The 
smaller fluctuations for large times reflect the intrinsic 
probabilistic nature of the model and are not related at 
all to the appearance of avalanches. This is clearly seen 
in Fig. H where we have plotted the cumulative number 



of aftershocks versus time instead of the differential plot 
of Fig. f§ As can be observed in Fig. [|, sudden accelera- 
tions appear in the first stages of rupture. This behavior 
closely resembles that previously reported in Sec. 2 for 
the eastern Pyrenees aftershock sequence (see Fig. [I]). 
The Omori's law maximum likelihood fit to the model 
results is also included in Fig. B. The estimated param- 
eter values are: p = 1.01 ± 0.06, K = 210 ± 29 and 
c = 1.3 ±0.2. 

In our model, the changes in aftershock rate are related 
to the readjustments of local stresses when events take 
place. With time, local concentrations of stress appear in 
the system and there is a high probability of finding a re- 
gion in which the load supported by the elements is close 
to the threshold value a — 1. That is, there is a large het- 
erogeneous stress state in which one failure will trigger 
an avalanche. During the evolution of the avalanche the 
local accumulation of load increases. This fact together 
with the high value of p, provokes the 5s corresponding 
to these steps of rupture to be considerably reduced. As 
a result, we observe the step-like change in the cumu- 
lative number of events (contrast Fig. 5 with Fig. 6). For 
large times the avalanches eventually disappear since in a 
non-conservative model the total load in the system sys- 
tematically decreases and hence it would be unlikely to 
accumulate stress in local regions as to surpass the value 
a= 1. 

Of interest is the further investigation of the acceler- 
ation events in order to get an additional insight about 
the observed aftershock sequences. One simple way to do 
that was explained in Sec. 2. It consists of decomposing 
the original series of aftershocks in leading events and 
cascades depending on whether a relaxation law is ac- 
complished or not. We have followed the same procedure 
with the synthetic data. The series of cascades obtained 
in such a way is shown in Fig. |[ Part (a) shows the series 
after removing all the cascades, i.e., leaving only leading 
aftershocks, and in Fig. |^b we plot the cascades, in which 
the first event of each cascade is a leading aftershock. The 
decomposition obtained from the model is indeed indis- 
tinguishable from that corresponding to the real series 
of events (Fig. ||). Two characteristics of the series of 
cascades are again relevant: one, the elapsed time be- 
tween successive leading events is larger than the preced- 
ing one in complete agreement with the supposition of a 
relaxation process, and second, the series of cascades can 
be well approximated by almost straight segments whose 
slopes decrease as time passes. This later characteris- 
tic could be used to quantify the observed jumps in the 
cumulative plot of aftershocks and to explain why they 
are present mainly in the first stages of rupture. The 
larger jumps are related to the occurrence of avalanches, 
which are caused by local accumulations of stress; so it is 
expected that when avalanches damp out due to dissipa- 
tion, changes in the rate of occurrence are more spaced 
in time and cascades consist of fewer events. Of course, 
there will be fluctuations about the power law trend even 
in the case where avalanches have deceased. Thus, we ex- 
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pect slope values gradually closer to zero as time tends 
to infinity. This is clearly appreciated in Fig. 0, where we 
have represented in a log-log plot the slopes of the cas- 
cades versus the occurrence time of the leading event that 
initiate each cascade. Plot (a) shows only the first part 
of the time sequence to facilitate comparison with Fig. ||. 
As for the observed series of aftershocks, the slopes fit a 
power law with an exponent of about v = 0.94 for the 
fist part of the time series and v = 1.08 for the long- 
tail end (plot (b)). Actually, the long-tail exponent v 
ranges between 1.00 and 1.08 depending on the conser- 
vation level 7r and the value of p. This appears to be 
a smooth dependence. Thus, the qualitative behavior 
is again captured. The discrepancy between the slopes 
(y = 0.7 for the Pyrenees sequence) is not surprising due 
to the simplicity of the model as compared with the in- 
herent complexity of the real phenomenon we want to 
simulate. The reason for this particular behavior, that 
is, why the slopes follow a power law and no other law, 
is unclear to us up to now. 

Future efforts will be devoted to the understanding of 
other dynamical characteristics of the model and their 
fine dependence on p and it by studying another com- 
plex series of aftershocks. We also plan to perform a 
detailed analysis of the spatial structure of the sequence 
of events coming from our model. 
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FIG. 1. Aftershock sequence of February 18, 1996, east- 
ern Pyrenees, (a) Complete series of aftershocks, shown as 
the accumulated number of events (left axis), together with 
their magnitude (right axis), (b) First 300 hours of the af- 
tershock sequence, as used in the comparison with the model 
results. The fit to the whole range and two independent fits 
to the 0-100 h and 140-300 h intervals has been superimposed 
(c) Separation of the aftershock sequence into leading after- 
shocks (filled circles) and cascade events (dots). See the text 
for details. 

FIG. 2. (a) Series formed by the leading aftershocks, after 
removal from the original sequence the cascades. The fit to 
the Omori law is much better than the original, and the p 
value (0.94) is also closer to worldwide aftershock p-values. 
(b) Cascades retrieved from the original first 300 hours of the 
aftershock sequence. Each cascade can be approximated by a 
straight line. 

FIG. 3. Slope of the cascades versus time on log-log scale. 
It can be clearly seen that it follows a power law s tx t~", 
with v ~ 0.7 

FIG. 4. Rate of aftershocks dN/dt as a function of dimen- 
sionless time for a dissipation of n — 0.7 and a Weibull index 
of p = 30. The spikes that decorate the general t -1 trend cor- 
respond to sudden accelerations in event rate (avalanches). 
The diagonal straight line has a slope of —1. This curve was 
obtained by numerically differentiating the curve in Fig. U 



FIG. 5. Accumulated number of aftershocks iV as a func- 
tion of dimensionless time for a dissipation of tt = 0.7 and a 
Weibull index of p — 30. Note the sudden increases in event 
rate (step-like jumps) superimposed to the general Omori-law 
trend. The continuous line is a maximum likelihood fit to the 
2-10 range and the corresponding p-value is p — 1.01 ± 0.06. 

FIG. 6. (a) Leading aftershock sequence for a simulation 
with 7r = 0.7 and p — 30. (b)Model cascades. The first event 
in each cascade is a leading aftershock. Note that the cas- 
cades can be also approximated by straight lines, as was the 
case with the cascades in the actual aftershock sequence, Fig. 
2. 

FIG. 7. (a) Slopes of the model cascades versus dimen- 
sionless time of the leading event that initiates each cascade 
on log-log scale. Only the first part of the model cascades is 
shown in this plot to facilitate comparison with Fig. 3. The 
local power-law exponent in this part is v = 0.94. (b) Log-log 
representation of the slopes versus time for the long-time tail 
of the simulation. As for the eastern Pyrenees series of after- 
shocks, the slopes fit very well a power law, in this case with 
an exponent of about 1.08. 
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